
clear all
set more off

use "merged.dta"

* table_1
label var percent_100 			"Final Course Grade (Percent Scale)"
label var temp_daily_avg 		"Temperature (Daily Avg. C)"
label var precip_daily_avg 		"Precipitation (Daily Avg. cm)"
label var relhum_daily_avg 		"Relative Humidity (Percent Scale)"
label var press_daily_avg 		"Pressure (kPa)"
label var winddir_cos_daily_avg	"East-West Wind Speed (Daily Avg. km/h)"
label var winddir_sin_daily_avg	"North-South Wind Speed (Daily Avg. km/h)"
label var pm25_days_above_10	"Number of Days in Past 90 with PM2.5$>$10"
label var pm25_days_above_12	"Number of Days in Past 90 with PM2.5$>$12"
label var pm25_semester			"90 Day Average PM2.5"
label var pm25_semester_med		"90 Day Median PM2.5"
label var co_daily_avg			"Daily Average CO (ppm)"
label var no2_daily_avg			"Daily Average NO2 (ppb)"
label var o3_daily_avg			"Daily Average O3 (ppb)"
label var so2_daily_avg			"Daily Average SO2 (ppb)"
label var inv_before_exam 		"Inversion Before Exam"
label var stem_course 			"STEM Course"
label var female_student 		"Female"
label var newbuilding 			"New Building (Built After 2000)"
label var leed 					"LEED Certification (Any)"
label var leed_gold 			"LEED Gold Certification"
label var upperfloor 			"Upper Floor (Above Ground Level)"

eststo clear
estpost tabstat ///
percent_100 ///
pm25_pref ///
pm25_semester ///
pm25_semester_med ///
temp_daily_avg ///
precip_daily_avg ///
relhum_daily_avg ///
press_daily_avg ///
winddir_cos_daily_avg ///
winddir_sin_daily_avg ///
pm25_days_above_10 ///
pm25_days_above_12 ///
co_daily_avg ///
no2_daily_avg ///
o3_daily_avg ///
so2_daily_avg ///
female_student ///
inv_before_exam ///
stem_course ///
newbuilding ///
leed ///
leed_gold ///
upperfloor ///
, listwise statistics(mean sd) columns(statistics)
esttab, cells("mean sd") nomtitle nonumber label
esttab using table_1.tex, replace cells("mean sd") nomtitle nonumber label

**********************************************************************************************************

* table_2_a
local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: reghdfe percent pm25_pref 				i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot) 
eststo: reghdfe percent pm25_pref `controls1'  i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot) 
eststo: reghdfe percent pm25_pref `controls2'  i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot) 

esttab using table_2_a.tex,  ///
keep(*pm25*) ///
mgroups("Preferred Specification",pattern(1 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_2_b
local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: reghdfe percent_no_scale pm25_pref 			i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot) 
eststo: reghdfe percent_no_scale pm25_pref `controls1' i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot) 
eststo: reghdfe percent_no_scale pm25_pref `controls2' i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot) 

esttab using table_2_b.tex,  ///
keep(*pm25*) ///
mgroups("No Scaling",pattern(1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_2_c
local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: reghdfe percent pm25_pref 				i.slot i.dow , absorb(strm) vce(clust stud_id exam_slot) 
eststo: reghdfe percent pm25_pref `controls1'  i.slot i.dow , absorb(strm) vce(clust stud_id exam_slot) 
eststo: reghdfe percent pm25_pref `controls2'  i.slot i.dow , absorb(strm) vce(clust stud_id exam_slot) 

esttab using table_2_c.tex,  ///
keep(*pm25*) ///
mgroups("No Student Fixed Effects",pattern(1 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table 3
* figure 4

capture drop pm25_bin
gen pm25_bin = .
replace pm25_bin = 0 if pm25_pref!=. & (pm25_pref >= 0 & pm25_pref<5)
replace pm25_bin = 1 if pm25_pref!=. & (pm25_pref >= 5 & pm25_pref<10)
replace pm25_bin = 2 if pm25_pref!=. & (pm25_pref >= 10 & pm25_pref<15)
replace pm25_bin = 3 if pm25_pref!=. & (pm25_pref >= 15 & pm25_pref<20)
replace pm25_bin = 4 if pm25_pref!=. & (pm25_pref >= 20)

capture label define pm25_bin_label 0 "0-" 1 "5-" 2 "10-" 3 "15-" 4 "20+"
capture label values pm25_bin pm25_bin_label
capture label define zero_label 0 "0-" 1 "5-" 2 "10-" 3 "15-" 4 "20+"
capture label values zero zero_label

local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

capture gen zero = 0

eststo clear

reghdfe percent zero i.pm25_bin i.slot i.dow i.strm , absorb(stud_id) vce(clust stud_id exam_slot) 
estimates store model1

reghdfe percent zero i.pm25_bin i.slot i.dow i.strm `controls1', absorb(stud_id) vce(clust stud_id exam_slot) 
estimates store model2

reghdfe percent zero i.pm25_bin i.slot i.dow i.strm `controls2', absorb(stud_id) vce(clust stud_id exam_slot) 
estimates store model3

esttab model1 model2 model3 using table_3.tex,  ///
keep(*pm25_bin*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
compress se(3) b(3) replace nogaps noomitted nomtitles ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

coefplot  (model1, label("No Controls"))(model3, label("With Weather Controls and Interactions")  msymbol(S) ciopts(lpattern(line)) ), keep(*zero* *pm25*) vertical omit ///
graphregion(color(white)) bgcolor(white) ///
xtitle("Final Exam PM2.5") xscale( titlegap(3)) xline() xlabel(1 "0-" 2 "5-" 3 "10-" 4 "15-" 5 "20+") ///
ytitle("Final Grade (Standard Deviations)") yscale(titlegap(3)) yline(0)  ylabel(, nogrid)
graph export "figure_4.pdf", as(pdf) replace


**********************************************************************************************************


*** table_4_a

local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear

	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_14_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 14	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day+pm25_lag_9_day+pm25_lag_10_day+pm25_lag_11_day+pm25_lag_12_day+pm25_lag_13_day+pm25_lag_14_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_13_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 13	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day+pm25_lag_9_day+pm25_lag_10_day+pm25_lag_11_day+pm25_lag_12_day+pm25_lag_13_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_12_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 12	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day+pm25_lag_9_day+pm25_lag_10_day+pm25_lag_11_day+pm25_lag_12_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_11_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 11	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day+pm25_lag_9_day+pm25_lag_10_day+pm25_lag_11_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_10_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 10	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day+pm25_lag_9_day+pm25_lag_10_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_9_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 9	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day+pm25_lag_9_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_8_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 8	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day+pm25_lag_8_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_7_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 7	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day+pm25_lag_7_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_6_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 6	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day+pm25_lag_6_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_5_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 5	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day+pm25_lag_5_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_4_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 4	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day+pm25_lag_4_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_3_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 3	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day+pm25_lag_3_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_2_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 2	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day+pm25_lag_2_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
	eststo : reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_1_day i.slot i.strm i.dow `controls2' , absorb(stud_id) vce(clust stud_id exam_slot) 
	estadd scalar maxlag = 1	
	
	lincom pm25_pref
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 

	lincom pm25_pref+pm25_lag_1_day
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	

esttab using table_4_a.tex ,  ///
drop(*) ///
mgroups("Unrestricted Finite Linear Lag",pattern(1 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label ///
stats(maxlag N  cum_b cum_s cum_t imp_b imp_s imp_t, ///
fmt(%9.0gc %9.0gc %9.2f ) ///
labels("Lags Included" "Observations" "Cumulative Effect" "S.E." "t" "Impact Effect" "S.E." "t" ) ///
layout("@" "@" "@" "(@)" "@" "@" "(@)" "@")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Semester FE = *strm*" ///
"Day-of-Week FE = *dow*" ///
"Time-of-Day FE = *slot*" ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#*" ///
, labels("Y" " ")) ///
nonotes	
beep


*** table_4_b

local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear

capture drop pdl_ml_*
forval s = 0(1)2 {
	gen pdl_ml_1_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day  
	gen pdl_ml_2_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day  	
	gen pdl_ml_3_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day  
	gen pdl_ml_4_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day  
	gen pdl_ml_5_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day  
	gen pdl_ml_6_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day  
	gen pdl_ml_7_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day	
	gen pdl_ml_8_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day  
gen pdl_ml_9_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day  
gen pdl_ml_10_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day  
gen pdl_ml_11_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day  
gen pdl_ml_12_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day  
gen pdl_ml_13_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day + 13^`s'*pm25_lag_13_day  
gen pdl_ml_14_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day + 13^`s'*pm25_lag_13_day + 14^`s'*pm25_lag_14_day 
}

forval i = 14(-1)1{
	capture drop mydeg0
	capture drop mydeg1
	capture drop mydeg2
	gen mydeg0 = pdl_ml_`i'_d_0
	gen mydeg1 = pdl_ml_`i'_d_1
	gen mydeg2 = pdl_ml_`i'_d_2

	reghdfe percent mydeg0 mydeg1 mydeg2 `controls2' i.strm i.dow i.slot, absorb(stud_id) vce(clust stud_id exam_slot) 
	constraint 1 mydeg0 + (`i'+1) * mydeg1 + (`i'+1) * (`i'+1) * mydeg2 = 0
	eststo: linest, c(1)  
	
	lincom mydeg0 + 0 * mydeg1 + 0* mydeg2
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 
	
	local j = (`i'*(`i'+1))/2
	
	local k = (`i'*(`i'+1)*(2*`i'+1))/6
	
	di `j'
	di `k'
	
	lincom (`i'+1)*mydeg0+(`j')*mydeg1+(`k')*mydeg2
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	

}

esttab using table_4_b.tex ,  ///
drop(*) ///
mgroups("Restricted Finite Distributed Quadratic Lag",pattern(1 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label ///
stats(maxlag N  cum_b cum_s cum_t imp_b imp_s imp_t, ///
fmt(%9.0gc %9.0gc    %9.2f  ) ///
labels("Lags Included" "Observations" "Cumulative Effect" "S.E." "t" "Impact Effect" "S.E." "t" ) ///
layout("@" "@" "@" "(@)" "@" "@" "(@)" "@")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
nonotes	



* table_4_c

local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear


capture drop pdl_ml_*
forval s = 0(1)3 {
	 gen pdl_ml_1_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day  
	 gen pdl_ml_2_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day  
	 gen pdl_ml_3_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day  	
	 gen pdl_ml_4_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day  
	 gen pdl_ml_5_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day  
	 gen pdl_ml_6_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day  
	 gen pdl_ml_7_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day
	 gen pdl_ml_8_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day  
gen pdl_ml_9_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day  
gen pdl_ml_10_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day  
gen pdl_ml_11_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day  
gen pdl_ml_12_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day  
gen pdl_ml_13_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day + 13^`s'*pm25_lag_13_day  
gen pdl_ml_14_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day + 13^`s'*pm25_lag_13_day + 14^`s'*pm25_lag_14_day 
}

forval i = 14(-1)1 {
	capture drop mydeg0
	capture drop mydeg1
	capture drop mydeg2
	capture drop mydeg3
	gen mydeg0 = pdl_ml_`i'_d_0
	gen mydeg1 = pdl_ml_`i'_d_1
	gen mydeg2 = pdl_ml_`i'_d_2
	gen mydeg3 = pdl_ml_`i'_d_3
	
	eststo: reghdfe percent mydeg0 mydeg1 mydeg2 mydeg3  `controls2' i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot)

	
	estadd scalar maxlag = `i'
	
	lincom mydeg0 + 0 * mydeg1 + 0* mydeg2
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 
	
	local j = (`i'*(`i'+1))/2

	local k = (`i'*(`i'+1)*(2*`i'+1))/6
	
	local l = (`i'^2)*((`i'+1)^2)/4
	
	lincom (`i'+1)*mydeg0+(`j')*mydeg1+(`k')*mydeg2+(`l')*mydeg3
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
}

esttab using table_4_c.tex ,  ///
drop(*) ///
mgroups("Restricted Finite Distributed Cubic Lag",pattern(1 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label ///
stats(maxlag N  cum_b cum_s cum_t imp_b imp_s imp_t, ///
fmt(%9.0gc %9.0gc    %9.2f  ) ///
labels("Lags Included" "Observations" "Cumulative Effect" "S.E." "t" "Impact Effect" "S.E." "t" ) ///
layout("@" "@" "@" "(@)" "@" "@" "(@)" "@")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
nonotes	
beep



* table_4_d

local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear


capture drop pdl_ml_*
forval s = 0(1)4 {
	 gen pdl_ml_1_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day  
	 gen pdl_ml_2_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day  
	 gen pdl_ml_3_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day  	
	 gen pdl_ml_4_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day  
	 gen pdl_ml_5_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day  
	 gen pdl_ml_6_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day  
	 gen pdl_ml_7_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day	
	 gen pdl_ml_8_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day  
gen pdl_ml_9_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day  
gen pdl_ml_10_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day  
gen pdl_ml_11_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day  
gen pdl_ml_12_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day  
gen pdl_ml_13_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day + 13^`s'*pm25_lag_13_day  
gen pdl_ml_14_d_`s' = 0^`s'*pm25_pref + 1^`s'*pm25_lag_1_day + 2^`s'*pm25_lag_2_day + 3^`s'*pm25_lag_3_day + 4^`s'*pm25_lag_4_day + 5^`s'*pm25_lag_5_day + 6^`s'*pm25_lag_6_day + 7^`s'*pm25_lag_7_day ///
	+ 8^`s'*pm25_lag_8_day + 9^`s'*pm25_lag_9_day + 10^`s'*pm25_lag_10_day + 11^`s'*pm25_lag_11_day + 12^`s'*pm25_lag_12_day + 13^`s'*pm25_lag_13_day + 14^`s'*pm25_lag_14_day 
}

forval i = 14(-1)1 {
	capture drop mydeg0
	capture drop mydeg1
	capture drop mydeg2
	capture drop mydeg3
	capture drop mydeg4
	gen mydeg0 = pdl_ml_`i'_d_0
	gen mydeg1 = pdl_ml_`i'_d_1
	gen mydeg2 = pdl_ml_`i'_d_2
	gen mydeg3 = pdl_ml_`i'_d_3
	gen mydeg4 = pdl_ml_`i'_d_4
	
	eststo: reghdfe percent mydeg0 mydeg1 mydeg2 mydeg3 mydeg4  `controls2' i.slot i.dow i.strm, absorb(stud_id) vce(clust stud_id exam_slot)
	
	estadd scalar maxlag = `i'
	
	lincom mydeg0 + 0 * mydeg1 + 0* mydeg2 + 0*mydeg3+0*mydeg4
	estadd scalar imp_b = r(estimate)
	estadd scalar imp_s = r(se)
	estadd scalar imp_t = r(estimate)/r(se) 
	
	local j = (`i'*(`i'+1))/2

	local k = (`i'*(`i'+1)*(2*`i'+1))/6
	
	local l = (`i'^2)*((`i'+1)^2)/4
	
	local m = ( (`i')*(`i'+1)*(2*`i'+1)*(3*(`i'^2)+3*(`i')-1)/30)         
	
	lincom (`i'+1)*mydeg0+(`j')*mydeg1+(`k')*mydeg2+(`l')*mydeg3+(`m')*mydeg4
	estadd scalar cum_b = r(estimate)
	estadd scalar cum_s = r(se)
	estadd scalar cum_t = r(estimate)/r(se)	
	
}

esttab using table_4_d.tex ,  ///
drop(*) ///
mgroups("Restricted Finite Distributed Quartic Lag",pattern(1 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label ///
stats(maxlag N  cum_b cum_s cum_t imp_b imp_s imp_t, ///
fmt(%9.0gc %9.0gc    %9.2f  ) ///
labels("Lags Included" "Observations" "Cumulative Effect" "S.E." "t" "Impact Effect" "S.E." "t" ) ///
layout("@" "@" "@" "(@)" "@" "@" "(@)" "@")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
nonotes	
beep



* table_5_a

local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: ivreghdfe pm25_pref inv_before_exam 			 if percent!=., absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  
eststo: ivreghdfe pm25_pref inv_before_exam `controls1'  if percent!=., absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  
eststo: ivreghdfe pm25_pref inv_before_exam `controls2'  if percent!=., absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  

esttab, keep(*inv_before_exam*)

esttab using table_5_a.tex,  ///
keep(*inv_before_exam*) ///
mgroups("First Stage",pattern(1 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_5_b

local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: ivreghdfe percent (pm25_pref = inv_before_exam)				, absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  
	estadd scalar CraggDonaldF = e(cdf)
	estadd scalar KleibergenPaapF = e(rkf)
eststo: ivreghdfe percent (pm25_pref = inv_before_exam) `controls1'	, absorb(slot dow strm stud_id) cluster(stud_id exam_slot) 
	estadd scalar CraggDonaldF = e(cdf)
	estadd scalar KleibergenPaapF = e(rkf)
eststo: ivreghdfe percent (pm25_pref = inv_before_exam) `controls2'	, absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  
	estadd scalar CraggDonaldF = e(cdf)
	estadd scalar KleibergenPaapF = e(rkf)
	
esttab, keep(*pm25_pref*)

esttab using table_5_b.tex,  ///
keep(*pm25_pref*) ///
mgroups("Second Stage",pattern(1 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N CraggDonaldF KleibergenPaapF, fmt(%9.0fc) labels("Observations" "CraggDonaldF" "KleibergenPaapF")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_5_c

local controls1 "temp_daily_avg precip_daily_avg relhum_daily_avg press_daily_avg winddir_cos_daily_avg winddir_sin_daily_avg"
local controls2 "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: ivreghdfe percent inv_before_exam			 , absorb(slot dow strm stud_id) cluster(stud_id exam_slot) 
eststo: ivreghdfe percent inv_before_exam `controls1', absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  
eststo: ivreghdfe percent inv_before_exam `controls2', absorb(slot dow strm stud_id) cluster(stud_id exam_slot)  

esttab, keep(*inv_before_exam*)

esttab using table_5_c.tex,  ///
keep(*inv_before_exam*) ///
mgroups("Reduced Form",pattern(1 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span ) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_6

eststo clear
local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"
			
eststo: reghdfe percent pm25_pref `controls' 								i.slot i.dow i.strm , absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref c.pm25_pref#1.female `controls' 			i.slot i.dow i.strm , absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref c.pm25_pref##1.stem_course `controls' 	i.slot i.dow i.strm , absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref c.pm25_pref##1.stem_course `controls' 	i.slot i.dow i.strm if stem_student==1, absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref c.pm25_pref##1.stem_course `controls' 	i.slot i.dow i.strm if stem_student==0, absorb(stud_id) vce(clust stud_id exam_slot)

esttab, keep(*pm25*)

esttab using table_6.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles("All" "All" "All" "STEM" "Non-STEM")    ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes	

* table_7

eststo clear
local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo: reghdfe percent pm25_pref pm25_semester  		`controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_semester_med   	`controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_days_above_10   `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_days_above_12   `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)

esttab, keep(*pm25*)

esttab using table_7.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
nobase mtitles() nomtitles  ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_8
		
eststo clear
local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo: reghdfe percent c.pm25_pref  				`controls' 								i.strm i.dow i.slot , absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.upperfloor 	`controls' 								i.strm i.dow i.slot , absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.upperfloor  `controls' 								i.strm i.dow i.slot if newbuilding==0, absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.upperfloor  `controls' 								i.strm i.dow i.slot if newbuilding==1, absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.newbuilding `controls' 								i.strm i.dow i.slot , absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.newbuilding c.pm25_pref##1.leed  		`controls' 	i.strm i.dow i.slot , absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.newbuilding c.pm25_pref##1.leed_gold  	`controls' 	i.strm i.dow i.slot , absorb(crse_id) vce(clust crse_id)
eststo: reghdfe percent c.pm25_pref##1.newbuilding c.pm25_pref##1.leed_gold  	`controls' 	i.strm i.dow i.slot if faculty!="SOC"  , absorb(crse_id) vce(clust crse_id)

esttab, keep(*pm25*) se

esttab using table_8.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle("\shortstack{Pref.\\Spec.}" "\shortstack{All \\ Bldg.}" "\shortstack{Old \\ Bldg.}" "\shortstack{New \\ Bldg.}" "\shortstack{All \\ Bldg.}" "\shortstack{Leed \\ Bldg.}" "\shortstack{Leed Gold \\ Bldg.}") ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes


* table_9

eststo clear

local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

capture drop pm25_custom
gen pm25_custom = pm25_pref
label var pm25_custom "PM2.5"

eststo: reghdfe percent pm25_custom  `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)

capture drop pm25_custom
gen pm25_custom = pm25_pref
replace pm25_custom = pm25_9 if hour==9
replace pm25_custom = pm25_14 if hour==14
replace pm25_custom = pm25_19 if hour==19
label var pm25_custom "PM2.5"

eststo: reghdfe percent pm25_custom  `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)

capture drop pm25_custom
gen pm25_custom = pm25_pref
replace pm25_custom = (pm25_7+pm25_8+pm25_9)/3 if hour==9
replace pm25_custom = (pm25_12+pm25_13+pm25_14)/3 if hour==14
replace pm25_custom = (pm25_17+pm25_18+pm25_19)/3 if hour==19
label var pm25_custom "PM2.5"

eststo: reghdfe percent pm25_custom  `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)

capture drop pm25_custom
gen pm25_custom = pm25_pref
replace pm25_custom = (pm25_9+pm25_10+pm25_11)/3 if hour==9
replace pm25_custom = (pm25_14+pm25_15+pm25_16)/3 if hour==14
replace pm25_custom = (pm25_19+pm25_20+pm25_21)/3 if hour==19
label var pm25_custom "PM2.5"

eststo: reghdfe percent pm25_custom  `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)

capture drop pm25_custom
gen pm25_custom = pm25_pref
replace pm25_custom = (pm25_1+pm25_2+pm25_3+pm25_4+pm25_5+pm25_6+pm25_7)/7
label var pm25_custom "PM2.5"

eststo: reghdfe percent pm25_custom  `controls' i.strm i.dow i.slot , absorb(stud_id) vce(clust stud_id exam_slot)

capture drop pm25_custom
gen pm25_custom = pm25_pref
replace pm25_custom = alt_pm25_9 if hour==9
replace pm25_custom = alt_pm25_14 if hour==14
replace pm25_custom = alt_pm25_19 if hour==19
label var pm25_custom "PM2.5"

eststo: reghdfe percent pm25_custom  `controls' i.strm i.dow i.slot if year>=2013, absorb(stud_id) vce(clust stud_id exam_slot)

capture drop pm25_custom
gen pm25_custom = pm25_pref
label var pm25_custom "PM2.5"

esttab using table_9.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle("\shortstack{Pref}" "\shortstack{T}" "\shortstack{(T-3,T)}" "\shortstack{(T,T+3)}" "\shortstack{1h-7h}" "\shortstack{Alt. Mtr.}" "\shortstack{2013 On}") ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes



* table_10

local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo clear
eststo: reghdfe percent pm25_pref 				`controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref co_daily_avg 	`controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref no2_daily_avg `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref o3_daily_avg 	`controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref so2_daily_avg `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref co_daily_avg no2_daily_avg o3_daily_avg so2_daily_avg `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)
esttab, keep(*pm25*)

esttab using table_10.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle("\shortstack{Pref.\\Spec.}" "\shortstack{Incl.\\ CO}" "\shortstack{Incl.\\ NO2}" "\shortstack{Incl.\\ O3}" "\shortstack{Incl.\\ SO2}" "\shortstack{Incl.\\ Multi}") ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"CO = *co_daily_avg*" ///
"NO2 = *no2_daily_avg*" ///
"O3 = *o3_daily_avg*" ///
"SO2 = *so2_daily_avg*" ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_11
*ssc install winsor

eststo clear

local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

capture drop pm25_winsor
gen pm25_winsor = pm25_pref
label var pm25_winsor "PM2.5"

eststo: reghdfe percent pm25_winsor `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)

capture drop pm25_winsor
winsor pm25_pref, gen(pm25_winsor) p(0.10) highonly
label var pm25_winsor "PM2.5"

eststo: reghdfe percent pm25_winsor `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)

capture drop pm25_winsor
winsor pm25_pref, gen(pm25_winsor) p(0.10) lowonly
label var pm25_winsor "PM2.5"

eststo: reghdfe percent pm25_winsor `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)

capture drop pm25_winsor
winsor pm25_pref, gen(pm25_winsor) p(0.10) 
label var pm25_winsor "PM2.5"

eststo: reghdfe percent pm25_winsor `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)

esttab using table_11.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle("\shortstack{Pref.\\Spec.}" "\shortstack{Wins.\\ High}" "\shortstack{Wins.\\ Low}" "\shortstack{Wins.\\ Both}" ) ///
nobase ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_12

eststo clear

local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo: reghdfe percent pm25_pref `controls', absorb(slot dow strm stud_id) vce(clust stud_id exam_slot)

eststo: reghdfe percent pm25_pref `controls', absorb(slot dow strm stud_id) vce(clust stud_id)

eststo: reghdfe percent pm25_pref `controls', absorb(slot dow strm stud_id) vce(clust stud_id date)

esttab, keep(*pm25*)

esttab using table_clustering_ols.tex,  ///
keep(*pm25*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle("\shortstack{Clust.\\ Stud. x Slot}" "\shortstack{Clust. \\ Student}" "\shortstack{Clust.\\ Stud. x Date}" ) ///
compress se(5) b(5) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_a1_a

eststo clear
local absorb "slot dow strm stud_id"
local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_1_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_2_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_3_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_4_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_5_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_6_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lead_1_day-pm25_lead_7_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)

esttab, keep(*pm25*)

esttab using table_a1_a.tex,  ///
keep(*pm25_pref*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle( "\shortstack{1 Lead}"  "\shortstack{2 Leads}"  "\shortstack{3 Leads}"  "\shortstack{4 Leads}"  "\shortstack{5 Leads}"  "\shortstack{6 Leads}"  "\shortstack{7 Leads}") ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes

* table_a1_b

eststo clear
local absorb "slot dow strm stud_id"
local controls "c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.temp_daily_avg##c.precip_daily_avg c.temp_daily_avg##c.relhum_daily_avg c.temp_daily_avg##c.press_daily_avg c.temp_daily_avg##c.winddir_cos_daily_avg c.temp_daily_avg##c.winddir_sin_daily_avg c.precip_daily_avg##c.relhum_daily_avg c.precip_daily_avg##c.press_daily_avg c.precip_daily_avg##c.winddir_cos_daily_avg c.precip_daily_avg##c.winddir_sin_daily_avg c.relhum_daily_avg##c.press_daily_avg c.relhum_daily_avg##c.winddir_cos_daily_avg c.relhum_daily_avg##c.winddir_sin_daily_avg c.press_daily_avg##c.winddir_cos_daily_avg c.press_daily_avg##c.winddir_sin_daily_avg"

eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_1_day pm25_lead_1_day-pm25_lead_1_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_2_day pm25_lead_1_day-pm25_lead_2_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_3_day pm25_lead_1_day-pm25_lead_3_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_4_day pm25_lead_1_day-pm25_lead_4_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_5_day pm25_lead_1_day-pm25_lead_5_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_6_day pm25_lead_1_day-pm25_lead_6_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)
eststo: reghdfe percent pm25_pref pm25_lag_1_day-pm25_lag_7_day pm25_lead_1_day-pm25_lead_7_day `controls', absorb(`absorb') vce(clust stud_id exam_slot)

esttab, keep(*pm25*)

esttab using table_a1_b.tex,  ///
keep(*pm25_pref*) ///
star(* 0.10 ** 0.05 *** 0.01) label stats(N, fmt(%9.0fc) labels("Observations")) ///
mtitle( "\shortstack{1 L+L}"  "\shortstack{2 L+L}"  "\shortstack{3 L+L}"  "\shortstack{4 L+L}"  "\shortstack{5 L+L}"  "\shortstack{6 L+L}"  "\shortstack{7 L+L}") ///
compress se(3) b(3) replace nogaps noomitted ///
indicate( ///
"Weather Controls = temp_daily_avg" ///
"Weather Interactions = *#c.precip_daily_avg*" ///
, labels("Y" " ")) ///
nonotes


* figure_1

hist pm25_pref if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Exam Day PM{sub:2.5}") xscale(titlegap(2)) xline(10) percent width(5) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_1.eps", as(eps) replace

* figure_2

hist pm25_days_above_10 if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Semester Days Above 10 PM{sub:2.5}") xscale(titlegap(2)) xline() percent width(5) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_2_a.eps", as(eps) replace

hist pm25_semester if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Semester Average PM{sub:2.5}") xscale(titlegap(2)) xline() percent width(1) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_2_b.eps", as(eps) replace

* figure_3

capture drop pm25_rounded
gen pm25_rounded = round(pm25_pref,0.1)

capture drop mean_percent
bysort pm25_rounded : egen mean_percent = mean(percent_100)

capture drop pm25_rounded_flag
egen pm25_rounded_flag = tag(pm25_rounded)

capture drop mean_percent_weight
bysort pm25_rounded : egen mean_percent_weight = count(percent)

twoway 	(scatter mean_percent pm25_rounded if mean_percent_weight>100 & pm25_rounded_flag [aweight = mean_percent_weight], mcolor(gs2) mlw(vthin) msize(small) msymbol(circle_hollow) yscale(titlegap(*10)) xscale(r(0 22) titlegap(*10))  ///
		ytitle("Exam Performance (%)") xtitle("Exam Day PM{sub:2.5}") xlabel(0(2)22) ylabel(66(1)76, nogrid) legend(off) graphregion(color(white)) bgcolor(white)  )
graph export "figure_3.pdf", as(pdf) replace

* figure_4 is above with the binned table

* figure_5
capture drop unique_date
egen unique_date = tag(year month day hour) 

capture drop density1 xaxis1 ci1_1 ci1_2
kdens pm25_pref if percent!=. & inv_before_exam==1 & unique_date , ll(0) bw(5) gen(density1 xaxis1) ci(ci1_1 ci1_2) nograph

capture drop density2 xaxis2 ci2_1 ci2_2
kdens pm25_pref if percent!=. & inv_before_exam==0 & unique_date, ll(0) bw(5) gen(density2 xaxis2) ci(ci2_1 ci2_2) nograph

twoway ///
(line density1 xaxis1, lpattern(solid) lcolor(black)) ///
(rarea ci1_1 ci1_2 xaxis1 , color(blue%10) lwidth(none none)) ///
(line density2 xaxis2, lpattern(dot) lcolor(black)) ///
(rarea ci2_1 ci2_2 xaxis2 , color(red%10) lwidth(none none)) ///
, graphregion(color(white)) ylabel(, nogrid) ///
legend(order(1 "Temperature Inversion" 3 "Normal Gradient") rows(2) position(0) bplacement(neast)) ///
xscale(titlegap(*10)) yscale(titlegap(*10)) ///
title("", color(black)) xtitle("Exam Day PM{sub:2.5}") ytitle("Density") 
graph export "figure_5.pdf", replace

* figure_a1
hist co_daily_avg if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Exam Day CO (ppm)") xscale(titlegap(2)) xline() percent width(0.05) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_a1_a.eps", as(eps) replace

hist no2_daily_avg if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Exam Day NO{sub:x} (ppb)") xscale(titlegap(2)) xline() percent width(5) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_a1_b.eps", as(eps) replace

hist o3_daily_avg if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Exam Day O{sub:3} (ppb)") xscale(titlegap(2)) xline() percent width(5) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_a1_c.eps", as(eps) replace

hist so2_daily_avg if percent!=., title() ///
ytitle("Percentage of Exams") ysc(titlegap(2)) ///
xtitle("Exam Day SO{sub:2} (ppb)") xscale(titlegap(2)) xline() percent width(0.5) start(0) lwidth(none) color(gs5) graphregion(color(white)) bgcolor(white)  ylabel(, nogrid)  gap(25)
graph export "figure_a1_d.eps", as(eps) replace

beep

